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Abstract. I investigate the process of converting gas into 
stars vifithin the framework of a standard cosmological 
model. By examining the set of objects grown in a com- 
bined N-body plus smoothed particle hydrodynamics sim- 
ulation with those obtained in similar models where some 
of the cold, dense gas was replaced by coUisionless "star" 
particles I show that it is possible to make this substi- 
tution without affecting the subsequent gas cooling rate. 
With even the most basic star forming criteria the masses 
of isolated objects are nearly identical to the mass of cold, 
dense gas found within the same objects in a non-star 
forming run. 

No evidence is found to support the contention that 
converting gas into stars might affect the amount of cold 
gas obtained in a simulation by retarding the cooling rate 
within those objects where stars have already formed. In 
practice, because cold gas can be reheated by shocks but 
stars remain as such whatever happens the masses of the 
largest objects found in the star forming runs are generally 
higher than those in the standard run. 

Finally, I demonstrate that an excellent match to the 
observed star formation rate can be achieved with even a 
very basic star formation prescription. 

Key words: methods; numerical - galaxies; formation; 
kinematics and dynamics - cosmology; theory - hydrody- 
namical simulation 



1. Introduction 

CoUisionless cosmological simulations have for several 
years been expanding our knowledge of how the large scale 
structure in the Universe formed. The difficulty lies in 
comparing the dark matter of the models with the observ- 
able Universe, where we see starlight from galaxies and X- 
rays from hot gas. Gas is a dissipational medium; tightly 
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bound structures form which survive the violent evolution 
of larger objects. When galaxies fall together to form clus- 
ters their dark matter halos become stripped and merge 
together, leaving behind several distinct, bright galaxies 
within a single dark matter halo. CoUisionless simulations 
lack the dissipational processes that occurred as visible 
galaxies cooled and formed stars which makes the process 
of comparing the results to the observations difficult. 

Recently several codes have been written that are 
capable of including a dissipating gaseous component 
(Katz, Hcrnquist & Weinberg 1992, Navarro & White 
1993, Evrard, Summers & Davis 1994, Tsai, Katz & 
Bertschinger 1994, Couchman, Thomas & Pearce 1995, 
Cen & Ostriker 1996, Frenk et al. 1996, Steinmetz 1996). 
Most of these groups employ smoothed particle hydrody- 
namics (SPH) to follow the gas (Monaghan 1992). This 
allows non-adiabatic heating through shocks and cooling 
via radiation to be included in the cosmological models 
and a direct comparison can be attempted between the 
regions of cold, dense gas and the observed distribution of 
galaxies. 

Here I examine the formation of objects within a cos- 
mological volume. Typically a Milky Way like galaxy will 
contain 100 or so gas particles with the mass resolution I 
employ. With this number of particles it is not possible to 
recover the internal dynamics and structure of the objects 
but reasonable distribution and multiplicity functions can 
be obtained. 

For a variety of reasons, not least because it saves com- 
putational effort, it is useful to convert the gas within cold, 
dense clumps (or "galaxies") into "star" particles, which 
from that point on behave like coUisionless dark matter 
particles. This can be achieved in a variety of ways, with 
varying degrees of physical motivation. Here I will attempt 
to show that perhaps the simplest possible choice of as- 
sumptions leads to a stable result and so more complicated 
schemes are perhaps unnecessary. Throughout this paper 
the term "star" refers to a particle that was once subject to 
gas forces but has since been converted into a coUisionless 
particle of the same mass. I do not consider any additional 
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physical processes which might affect a stellar population 

such as supernova feedback or metal enrichment as these 
will be discussed in detail elsewhere. 

Within this paper I show that one successful way of 
incorporating star formation within a cosmological simu- 
lation without affecting the masses of the final objects is 
to use a combined density and temperature cut. This al- 
leviates the worry that the action of converting a gaseous 
particle into a collisionless one might retard subsequent 
cooling by lowering the local gas density. In practice this 
is not the case and apart from dynamical considerations 
we are free to equate the mass of stars found with the 
mass of cold, dense gas that would have been obtained in 
a non-starforming simulation. 

I then compare the star formation rate obtained with 
my preferred method against that observed by Madau 
(1996). The good agreement between the observations and 
the simulations demonstrates that reproducing the ob- 
served star formation rate is in practice straightforward 
and is a natural consequence of hierarchical clustering. 

2. Simulations 

I have run my simulations using a parallel adaptive 
particle-particle, particle-mesh plus smoothed particle hy- 
drodynamics code (Pearcc & Couchman 1997) imple- 
mented in CRAFT, a directive based parallel Fortran de- 
veloped by Cray. It can be run in parallel on a Cray T3D or 
serially on a single processor workstation and is essentially 
identical in operation to the publicly released version of 
Hydra (Couchman, Pearce & Thomas 1996) which is de- 
scribed in detail by Couchman, Thomas & Pearce (1995). 
This work has been carried out as part of the programme 
of the Virgo Consortium, a group of astrophysicists inter- 
ested in large cosmological simulations. 

The three simulations presented here were of an flo = 
1, standard cold dark matter Universe with a boxsize 
of 25h~^ Mpc. I take h = 0.5 throughout this work, 
equivalent to a Hubble constant of 50kms~^ Mpc"^ . The 
baryon fraction, Clb was set from nucleosynthesis con- 
straints, rihh'^ = 0.015 (Copi, Schramm & Turner 1995) 
and I assume a constant gas metallicity of 0.5 times the 
solar value. Identical initial conditions were used in all 
cases, allowing a direct comparison to be made between 
the objects formed. Useful parameters for the runs are 
listed in table 1. 

The initial fluctuation amplitude was set by requir- 
ing that the model produced the same number density 
of rich clusters as observed today. To achieve this I take 
(Tg = 0.6, the present day linear rms fluctuation on a 
scale of 8/1-1 ^p^, ^^gj^^^ (j^i^ ^ p^^^y. ^qqq^ Yiana & Lid- 
die 1996). Each model began with 64^ dark matter par- 
ticles each of mass 3.1 x IO^^Mq and 64^ gas particles 
of mass 1.9 X 10^ M©, smaller than the critical mass de- 
rived by Steinmetz & White (1996) required to prevent 
2-body heating of the gaseous component by the heav- 



ier dark matter particles. I employ a comoving Plummer 

softening of 25h^^ kpc, which is typical for cosmological 
simulations but much larger than required to accurately 
simulate the dynamics of galaxies in dense environments. 
For a Plummer softening of 25h~^'kpc the effective soft- 
ening is around 50h~^ kpc whilst elliptical galaxies have 
scalelengths of around 4/1" ^ kpc. This dramatic increase in 
the cross section of the objects will artificially enhance the 
merger rate and make tidal stripping much more effective. 

With only 2 x 64^ particles converting the dense gas 
into collisionless particles does not resiilt in a large time 
saving. The total number of steps taken is not reduced 
by much because the dissipation that occurred as the gas 
cooled to form dense knots produces tightly bound clumps 
of stars which require small timesteps if they are not to be 
disrupted. A modest saving is seen in the CPU time re- 
quired per step as for stars expensive SPH calculations no 
longer need to be carried out. The full benefit of converting 
gas into stars only really becomes apparent within larger 
2 X 128^ particle simulations where shorter timesteps are 
observed (Pearce et al. 1997) if star formation is not em- 
ployed because the volume contains larger objects which 
form earlier. 

Before star formation can be introduced into a model 
an adequate treatment of gas cooling is required. In this 
paper I use the analytic form of Thomas & Couchman 
(1992). I have implemented the tabulated cooling func- 
tions of Sutherland & Dopita (1993) and no significant 
differences to the results presented here are obtained. In 
practice the results presented in this paper are not sensi- 
tive to the precise form of the cooling function and should 
hold for all such functions. 

The gas ends up in three phases; there is a hot phase, 
where gas sits in the potential wells formed by the larger 
objects, and two cold phases, one formed by the collapsed 
objects and the other formed by gas that has yet to col- 
lapse which occupies the void regions. Table 1 lists the 
percentage of the gas which resides in each of these phases 
at a redshift of zero. These values are in agreement with 
those obtained by Evrard, Summers & Davis (1994) at a 
redshift of 1 (the end of their simulation), where I obtain 
values of 24 percent hot gas and 72 percent cold gas in 
voids in all three simulations, the remaining 4 percent is 
in the form of cold, collapsed gas or stars. 

2.1. Star formation details 

Stars are formed within the models in two slightly differ- 
ent ways. In both cases once a gas particle has satisfied 
the relevant criteria its type is simply changed from "gas" 
to "star" and from then on the particle behaves like a col- 
lisionless dark matter particle of the same mass, only ex- 
periencing the force of gravity. This "instantaneous" con- 
version of a gas particle into a star is deliberately crude, 
reducing the number of free parameters to a minimum 
and is similar to that employed by Summers (1993). I do 
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R,uri 


AT 


AT 


/Og 


/Oh 


/Oc 


Steps 


cpu 


No stars 


31369 


315 


97 


46 


43 


1897 


20.8 


P + T 


29740 


327 


86 


45 


44 


1732 


16.6 


p + T + % 


29033 


320 


87 


45 


44 


1822 


18.1 



Table 1. Parameters for the 3 runs. These are the run without 
any star formation, one employing a density plus a temperature 
cut, and finally a density, temperature and small probability 
cut. The columns list the number of star particles at the end 
(or for the no star case the immber of cold, dense particles 
which would have formed stars if a density plus a temperature 
cut were employed on the final step), the number of groups 
found by the group finder at the end of the simulation, the 
percentage of the stars that arc within objects defined by the 
group finder, the percentage of the gas that is above 10'' A", the 
percentage of gas that is below 10* if and at low density, the 
number of system steps required and the total CPU time for 
the entire run in hours on 64 T3D processors. 



not consider schemes which involve spawning additional 
particles {e.g. Katz 1992, Navarro & White 1992), redis- 
tributing mass or dual identity "schizophrenic" particles 
(Mihos & Hernquist 1994). All these schemes may increase 
the tunability of the star forming process but as I will later 
show they arc not required at least in a cosmological sit- 
uation where the poor mass resolution makes any model 
of star formation very crude. 

In the first model a star is formed if a gas particle 
has a density exceeding 550 times the mean density in 
the box and has a temperature below 10''' if, converting 
11 percent of the gas into stars by the end of the sim- 
ulation {Vl^, = nb X 0.11 = 0.006). This is close to the 
observed stellar mass fraction, fi^ = 0.004 (Peebles 1993) 
despite the poor mass resolution of these simulations. This 
is another demonstration of the well known cooling catas- 
trophe (White & Frenk 1991). The second method is the 
same as the first except that a gas particle only has a 5% 
chance of becoming a star on any one step. An overden- 
sity of 550 is quite low compared to the threshold used 
by other authors {e.g. Navarro & White 1992) but I am 
forced to adopt this value because of my poor mass res- 
olution and large softening which causes the maximum 
reliable overdensity obtained within the simulation to be 
around 2000. In practice cold gas in this regime has little 
pressure support and so rapidly increases in density. 

Within a cosmological simulation it seems non-sensical 
to adopt a physical density criterion as the basis for a star 
formation algorithm. The star forming regions are charac- 
terised as collapsed clumps of cold, dense particles where 
the important factor is not the physical density but the 
overdensity. At a high enough redshift the entire Universe 
will be at a density above whatever physical threshold is 
employed, forcing the implementation of an overdensity 
constraint to prevent star formation at early times. If the 
star forming region itself could be resolved within a sim- 
ulation then a physical density threshold might be useful, 



but when an entire galaxy is only barely resolved such a 
threshold seems at best problematic. 

2.2. Identifying groups 

Within a simulation volume it is useful to be able to re- 
liably identify a catalogue of all the bound objects. In 
practice this is a tricky procedure because some degree of 
merging and disruption of the object set will be taking 
place at any given time and a binary merger can dramat- 
ically affect the position of any given object within the 
catalogue. The presence of a cooling gaseous component 
reduces the problem of identifying groups because both 
density and temperature information is available. The ob- 
jects being sought are characterised as clumps of cold, 
dense particles. In principle both density and tempera- 
ture information could be calculated for the dark matter 
component but the lack of dissipation and large scatter 
reduces the usefulness of the local velocity dispersion al- 
though others have used the local dark matter density as 
an aid to defining groups (Gelb & Bertschinger 1994). 

To define groups I use the procedure of Thomas & 
Couchman (1992). By extracting all the gaseous particles 
which are simultaneously both in a dense region and have 
a temperature below lO^K a reduced set of positions is 
obtained which can then be passed to a friends-of-friends 
algorithm. When the maximum linking length is small the 
particles within each clump are linked together. At higher 
values the recovered object set remains static until neigh- 
bouring clumps begin to be linked. This allows us to use a 
maximum linking length of 60h~^ kpc, a value larger than 
normal because of the highly reduced set of points that 
are passed to the finder. In practice the maximum link- 
ing length can vary over a large range and essentially the 
same object set be recovered. 

Once star formation is implemented the stars them- 
selves provide an additional method for making the initial 
selection of those particles which might be in groups. For 
the runs which contain star particles the friends-of-friends 
finder uses the star positions to produce the group cata- 
logue with the same linking parameters as for the gas runs. 
With stars the isolated objects are still easy to recognise 
and define but within dense regions things are much more 
messy. Clumps of stars are often non-spherical due to tidal 
distortions and the whole cluster is permeated by a dif- 
fuse background light that originates from small systems 
that have been completely disrupted. The central part of 
this halo is linked with the mass of the central object but 
in practice few galaxies are linked together because the 
stellar systems remain centrally concentrated. In total 86 
percent of the stars end up in the object set produced 
by the group finder. For the gas only run a much higher 
value of 97 percent is achieved, simply demonstrating the 
messier nature of runs which include stars. Practically, the 
merger and halo effects could be alleviated by reducing the 
value of the linking length used to define the stellar clus- 
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ters but this would further reduce the percentage of stars 
lying within resolved objects. 

3. Results 

3.1. Multiplicity function 

The group finder described above produces a catalogue 
of objects ordered by mass. The effect on the multiplicity 
function of forming stars using my different prescriptions 
is shown in figure 1. 

The gas particles being considered for conversion into a 
star should always be both dense and cold, residing within 
collapsed objects. If only a density threshold were em- 
ployed it would be necessary to choose a value low enough 
to allow small, isolated clumps to form stars but also high 
enoiigh to prevent hot gas particles in the cluster halo 
being converted. The extra computational effort of exam- 
ining the temperature is small and removes the problem of 
dense, hot cluster halos being converted into stars. As an 
artifact of the SPH approach it is also possible for a hot 
halo particle to overestimate its local density if it passes 
close to the central cold gas clump. With just a density 
cut these particles may be erroneously flagged as stars. 

The most massive stellar objects formed within the 
starforming runs arc generally twice as massive as those 
in the standard run. The clustered regions where these 
large objects reside are not being overmerged by the group 
finder as near the largest object the group finder recovers 
13 objects in the standard run as opposed to 12 in the 
stellar run. One tiny object near the centre of the star 
cluster has been merged, raising the central object mass by 
less than one percent. This shows that the central object 
really does contain twice as many stars as the amount 
of cold gas in the non-starforming nm. Dense, cold gas 
provides little pressure support and so its removal should 
not affect the cooling rate of the halo, as demonstrated by 
the observation that the central temperature of the hot 
halo when stars are present is only 1.2 percent below that 
of the standard (no star) run and the percentage of gas in 
the hot phase is almost constant between the three runs. 

The reason for the more massive stellar objects is that 
small objects can be tidally disrupted as they fall into a 
cluster (Summers 1996, Navarro & Steinmctz 1997). When 
the stellar clumps become tidally disrupted the stars re- 
main stars and form a diffuse halo, some of which gets 
attributed to the central object. 15 percent of those par- 
ticles which are both cold and dense enough to have been 
called stars at a redshift of 0.5 do not remain cold and 
dense at the end of the simulation. These particles pre- 
dominantly lie in objects that are in the process of in- 
falling into the larger clusters within the simulation vol- 
ume and so the disruption process is heavily biased. This 
disruption should be seen as a numerical artifact due to 
the large softening length employed. With a smaller soft- 
ening length both the cold, dense clumps and the stellar 



associations would have been allowed to dissipate further 
as they formed and perhaps survive passage through the 
harsh cluster environment. 

3.2. Star formation history 

The redshift at which stars were created is shown in fig- 
ure 2. The first "stars", formed in the precursor objects 
of what become the central galaxies of the largest clus- 
ters, appear around a redshift of 4. This initial formation 
epoch should be treated with extreme caution because if 
the resolution were to be increased then these first stars 
would be seen even earlier. This is the well known "cooling 
catastrophe" where more and more gas cools as the res- 
olution is progressively increased (White & Frenk 1991). 
The formation time of the earliest objects also depends 
upon the size of the volume being simulated; bigger boxes 
have more room to fit in a high sigma peak of the initial 
fluctuation spectrum, providing a site where a large clus- 
ter will eventually form. If the resolution is fixed bigger 
boxes will form their first objects earlier. 

As can be seen the standard, non-star forming run al- 
ways has more particles which would have qualified as 
stars if the same density and temperature cut were em- 
ployed as for the star forming runs. This is despite the 
disruption of objects mentioned in the previous section. 
The reason for this discrepancy is that the star forming 
algorithm is discrete, but the SPH variables used to mea- 
sure the state of the gas arc averages over some number, 
Ngph, of particles (here Ngph = 32). For a single star to 
be formed there must be Ngph particles close to the star 
formation threshold and so a zero point offset has to be 
introduced if the object masses are to be compared. To 
compensate for this effect I add aNsph onto the mass of 
all the stellar groups, where a is a parameter between 
and 1. For simplicity I take a = 1 in what follows. 

In figure 3 I plot the star formation rate in 
MQ/yr /Adpc? versus redshift with the observed results 
from Madau (1996) overlayed. Clearly stars are being over- 
produced at late times (after a redshift of 1) but this may 
just be a consequence of the cosmological model employed 
(standard CDM has lots of evolution at late times). This 
demonstrates that any star formation formula that is ul- 
timately based upon the gas density and is normalised to 
produce the correct total mass of stars should be expected 
to fit the observed star formation rate very well. 

As can be seen from figures 1, 2 & 3 the runs with a 
percentage chance of forming a star on any one step pro- 
duce almost identical results to their counterparts. The 
initial motivation for introducing a percentage chance to 
delay star formation was to assist gas cooling by pro- 
viding a "seed" of cold gas which was not immediately 
converted into stars onto which more gas could accumu- 
late. The success of the model with a combined density 
and temperature cut demonstrates that this concern was 
unfounded and that delaying star formation via a ran- 
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Fig. 2. The star forming history 
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Fig. 3. The star formation rate in Mq /yr /Mpc^ versus redshift for the simulations with a density and temperature cut with 
and without a delay via a percentage chance. The large crosses show the observed points taken from Madau (1996) 



dom chance is an unnecessary complication that in prac- 
tice simply raises the density threshold at which stars are 
formed because during the delay introduced the clump 
collapses further. This has the effect of producing more 
tightly bound clumps which may then survive subsequent 
tidal disruption but in practice has little effect upon the 
simulation. This lack of effect occurs for the same reason 
that a zero point offset is required; at any one time Ngph 
particles arc close to the star formation threshold so the 
addition of a few extra particles above the threshold has 
little effect. 

3.3. Direct comparison 

In figure 4 I show a comparison between the masses of the 
objects found by the group finder in the standard non-star 
forming run and the masses of the same objects found in 
my preferred star forming run, where both density and 
temperature cuts were employed but without a percentage 
chance (which has little effect). 

The correspondence in mass between the objects in 
the two runs is very good with only a small scatter. At the 
high mass end the masses of the objects in the star forming 
run are in fact higher than the equivalent object in the run 
without stars. This is because some objects falling into the 
cluster halos are disrupted and dispersed. When stars are 
present these subsequently add to the mass of the central 
object but without stars the gas is reheated and forms 
part of the hot cluster halo. With a smaller gravitational 
softening these objects would be more tightly bound and 



harder to disrupt and so the object masses in the non-star 
forming run would have been higher. 

The few outliers in the upper left quadrant of the fig- 
ure are due to objects within the star forming simulation 
that have been linked together by the group finder. This 
is because with stars the objects often have a small "halo" 
because they cannot dissipate further once stars have been 
formed and subsequent mergers and tidal torques lead to 
messier shapes, forming bridges and spurs which some- 
times cause the group finder to link together separate ob- 
jects. The paucity of objects in this region of the plot 
clearly demonstrates that overmerging is not a problem 
with these simulations. 

4. Conclusions 

The main conclusion of this work is that replacing cold, 
dense gas with stars within a cosmological simulation can 
be done with relative ease and to some degree of accuracy 
without affecting the subsequent cooling rate of the gas. 
In principle all that is required is that the gas is dense 
enough but in practice a temperature threshold should 
also be employed to ensure that hot halo particles are not 
converted to stars erroneously. 

Introducing a percentage chance to delay star forma- 
tion raises the density threshold and has no affect on the 
masses of the objects obtained. This and more complicated 
schemes such as the gradual conversion of a gas particle 
into a star might be useful if additional physical processes 
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Fig. 4. Comparison of the masses of objects with and without star formation. The hne shows where identical masses would fall. 



such as feedback are considered or if a smoother star for- 
mation rate is desired but are not required in principle. 

It is a fundamental property of the SPH method that 
if one particle has a certain density and temperature Ngph 
particles should have similar properties. This leads to a 
zero point offset in the masses of the stellar objects when 
compared to the mass of cold, dense gas in a non-star 
forming run. 

Matching the star formation rate obtained within a 
simulation to that observed is a straightforward exercise 
with good fits being obtained with even the most basic 
star formation formulae. This is because the general trend 
of the observed values, a small rate at early times rising 
slowly to a peak around a redshift of 1 followed by a de- 
cline just mimics the effects of structure formation in any 
underlying cosmology. 

Although the introduction of star formation allows 
simulations to be carried out more rapidly the simulator 
should be aware of several additional processes it intro- 



duces. Isolated objects are well reproduced whether or not 
star formation is employed but within dense regions ad- 
ditional dynamical processes are taking place. Obviously 
the dynamics of a gaseous clump will be different from 
a collisionless object in an environment where ram pres- 
sure stripping and merging are important. Small cold gas 
clumps can be disrupted as the cluster forms, dissolving 
into the general cluster halo. The analogous stellar sys- 
tems will also be tidally disrupted but the stars themselves 
will survive and contribute to a diffuse cloud of stellar par- 
ticles that pervades the larger groups. Gas clumps merge 
quite easily whereas star clumps tend to produce messy 
associations because once a star has been formed no sub- 
sequent dissipation can take place and so reheating events 
can prove to be a problem. 

In summary, converting gas into stars within a cos- 
mological N-body simulation is a viable option and even 
a basic prescription works well for isolated objects. How- 
ever, in regions where tidal torques and merging are sig- 
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nificant the differences between coUisionless systems and This article was processed by the author using Springer- Verlag 

collisional, dissipating gaseous systems should be carefully I^T^X A&A style file L-AA version 3. 
considered and if possible a parallel non-star forming sim- 
ulation should be carried out. 
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